Impurity-induced step interactions: a kinetic Monte-Carlo study 



Jiirgen Vollmer, 1, 2, Q Jozsef Hegediis, 1 ' 3 Frank Grosse, 4 and Joachim Krug 5 

1 Fachbereich Physik, Philipps Universitdt, Renthof 6, 35032 Marburg, Germany 
2 Dept. Dynamics of Complex Fluids, Max Planck Institute for Dynamics 
and Self-Organization, Bunsenstr. 10, D-37073 Gottingen, Germany 
1 Department of Chemistry, University of Cambridge, Lensfield Road, CB2 1EW Cambridge, UK 

*Institut fur Physik, Humboldt- Universitdt zu Berlin, Newtonstr. 15, 124-89 Berlin, Germany 
" 'Institut fur Theoretische Physik, Universitdt zu Koln, Ziilpicher Str.77, 50937 Koln, Germany 

(Dated: February 2, 2008) 

A one-dimensional continuum description of growth on vicinal surfaces in the presence of immobile 
impurities predicts that the impurities can induce step bunching when they suppress the diffusion of 
adatoms on the surface. In the present communication we verify this prediction by kinetic Monte- 
Carlo simulations of a two-dimensional solid-on-solid model. We identify the conditions where quasi 
one-dimensional step flow is stable against island formation or step meandering, and analyse in detail 
the statistics of the impurity concentration profile. The sign and strength of the impurity-induced 
step interactions is determined by monitoring the motion of pairs of steps. Assemblies containing 
up to 20 steps turn out to be unstable towards the emission of single steps. This behavior is traced 
back to the small value of the effective, impurity-induced attachment asymmetry for adatoms. An 
analytic estimate for the critical number of steps needed to stabilize a bunch is derived and confirmed 
by simulations of a one-dimensional model. 

PACS numbers: 81.15.Aa, 68.55.-a, 68.55.Ln 
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I. INTRODUCTION 



Impurities and adsorbates affect the growth of crystals 
and thin films in a variety of ways. Small amounts of 
CO strongly enhance the nucleation density in homoepi- 
taxial growth of Pt on Pt(lll) and reverse the orienta- 
tion of the resulting triangular islands [l|, 0] ; a floating 
monolayer of a suitably chosen surfactant species induces 
layer-by-layer growth in many growth systems 0, H[ ; and 
the adsorption of antifreeze proteins on the growth sur- 
face prevents the formation of macroscopic ice crystals in 
the blood of fish living in polar waters [J] . 

On a vicinal surface growing by step propagation, im- 
purities are generally expected to slow down the steps by 
pinning. This can lead to the formation of step bunches 
[H HI, Q- A different mechanism for impurity- induced 
step bunching not related to step pinning was recently 
proposed in the context of SiC growth on Si(100) where 
C plays the role of a codeposited impurity [1, 1| (see also 
Due to the motion of the steps, at any given time 
different parts of a terrace have been exposed to the im- 
purity flux for different durations, which leads to a gra- 
dient in the impurity concentration directed towards the 
ascending step. The coupling of the impurity concentra- 
tion profile to the diffusion of the growth units on the 
terrace may destabilise the equidistant step train. Orig- 
inally [| a lower binding energy to the impurities was 
suggested and confirmed by kinetic Monte Carlo (KMC) 
simulations as a possible physical source for the experi- 
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mentally observed behavior. Another scenario was dis- 
cussed in Q : Impurities slow down the adatom diffusion 
without affecting adatom binding energies (random bar- 
riers [ll|) which also leads to instability. 

The linear stability analysis of was based on a one- 
dimensional model of straight steps. The impurity and 
adatom concentration fields were treated in a continuum 
approximation and assumed to take on their stationary 
profiles instantaneously on the time scale of step mo- 
tion. In the present communication we revisit the prob- 
lem within a fully microscopic KMC simulation, taking 
explicit account of non-stationarity and fluctuations. 

We identify a range of parameters where the assump- 
tion of a one-dimensional array of straight steps is appli- 
cable. We then consider systems of two, three and more 
steps with periodic boundary conditions in order to fol- 
low the loss of stability of the equidistant arrangement, 
and characterise the long-term evolution of the system. 
The basic stability properties predicted by the continuum 
theory are confirmed. Simulations with two steps show 
bound pairs and equidistant steps with only small fluctu- 
ations of the terrace width in the appropriate parameter 
regimes. However, we also find that the impurity-induced 
step interactions are unable to stabilise larger assemblies 
of steps. Simulations with 3, 4, 8, and 20 steps approach 
highly dynamic states with many closely adjacent pairs 
of steps that frequently exchange partners. This behavior 
can be explained within the framework of a deterministic 
step-dynamical model [121 ] : in spite of attractive interac- 
tion between the steps, bunches that contain less than a 
critical number of steps decay by step emission. 

The paper is organized as follows. In the next sec- 
tion the KMC model is introduced and the fundamental 
growth modes are identified as a function of the system 
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parameter 

J = ex P 

Jim = exp 
V = 3/$ 



description 

suppression of diffusion by bonds 

change of diffusion by impurities 

hopping rate / deposition rate 

fraction of impurities in deposited atoms 



TABLE I: Dimensionless parameters characterising the diffu- 
sion of adatoms and surface growth. 



parameters. Section IIIII contains a detailed analysis of 
the spatial distribution of impurities on the terraces, fo- 
cusing in particular on the fluctuations around the mean 
impurity concentration gradient. The dynamics of step 
pairs, triplets and bunches is described and analyzed in 
Sec lIVl and conclusions are given in SecfVl 



II. MICROSCOPIC GROWTH MODEL 

We consider an SOS system with a simple cubic lattice, 
where the surface has an extension Lj x Lj. There are 
periodic boundary conditions along i, and Lees-Edwards 
boundary conditions along j. The latter are also periodic 
on the surface, but they induce a change of height of N 
steps when transversing the system in j-direction. By 
this topological constraint one enforces the existence of 
N steps on the surface which are aligned parallel to i. 
Depending on the growth parameters these steps may 
be fairly straight, they may meander, or there may be 
additional islands on the surface. 

The kinetic Monte Carlo simulation is constructed to 
be close to the continuum description Q . In the following 
the technical details are described with special emphasis 
on the differences to previous simulations Q. The tran- 
sition rate of a thermally activated hopping process from 
site A to the neighboring site B is given in transition 
state theory by 
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The preexponential factor To is taken usually to be 10 13 
s _1 , but other values are possible also [HI]. The activa- 
tion energy 
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is given by the difference of transition energy Et and 
binding energy Eb- Et may depend on initial and final 
state, whereas Eb only depends on the initial state. From 
now on these dependencies are suppressed. In the case 
of a simple cubic lattice usually the binding energy is 
described by a next neighbor counting model [lij 



with n being the number of in-plane next neighbors. In 
the present work the impurities are taken to influence the 
transition energy Et only, 



Et = E t + E- n 



(4) 



with Et being the transition energy of free adatoms. Pre- 
vious simulations Q had only considered the influence 
on the binding energy Eb, which leads to a strong step 
bunching effect. Inserting Eq. ([3]) and ([4} into Eq. ((T|) 
results in the hopping rate of an adatom residing on an 
impurity 
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The factors J and J; m arc defined in Tab. HI Therefore, 
free adatoms on the surface perform a random walk with 
hopping rate g. Their hopping rate is modified by the 
existence of neighbors and impurities. 

In the simulation adatoms are added randomly to the 
surface with a flux of $ atoms per lattice site, and the 
ratio of diffusion and incoming flux will henceforth be 
characterized by v = (7/$. An adatom on the surface will 
perform on average v/(LiLj) free steps before another 
atom is added anywhere on the surface. 

Impurities disappear from the surface by burying them 
under normal atoms, and they are created by a flux onto 
the surface. A fraction of p of impinging atoms are im- 
purities. When they hit the surface they immediately 
exchange position with a "normal" atom in the surface. 
Subsequently, there is a new impurity in the surface, and 
an additional adatom diffusing on the surface. 

Altogether the dynamics of surface growth is hence 
characterized by four dimensionless parameters summa- 
rized in Table HI Kinetic Monte Carlo simulation of this 
model (c/ [H} for details on the algorithm) show that 
the model is capable of reproducing the crossover from 
step flow at large where adatoms mostly attach to step 
edges, to island nucleation for smaller v, where adatoms 
merge into islands before reaching the step edges. More- 
over, the simulations also show the formation of step pairs 
for appropriately chosen binding strength J lm and den- 
sity of impurities p (see Sec lIVp . The different growth 
regimes for a system with N — 2 steps are summarized 
in Fig. [1] Fig. [2] shows snapshots of the time evolution 
of the surface height and the distribution of impurities 
for a system of two equidistant steps which evolves into 
a steady state where the two steps form a pair. 



III. DISTRIBUTION OF IMPURITIES 

Before further discussing the numerical results it is il- 
luminating to calculate the distribution of impurities on 
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FIG. 1: Summary of findings for the dominant growth mech- 
anism in a system with two steps as a function of Jj m and v. 
The other system parameters are fixed to the values Li — 25; 
Lj = 100; J = 40; p = 0.1. 



the terraces. To this end we consider two steps which 
are roughly aligned in parallel such that they enclose a 
terrace of width it). In a system of two parallclly aligned 
steps the width w of the terrace does not change in a 
steady state where the steps form a bound pair, such as in 
the lower-most panels of Fig. [2] A small surface element 
is created in this situation when the edge of the step is 
formed at that position and height, and it is buried after 
deposition of two ML, when the two steps have reached 
the same position again. In order to study the distribu- 
tion of impurities we use a comoving coordinate frame 
where i denotes the position in lateral direction, and d 
the distance from the position of the steps measured in 
the direction of growth. The latter distance will be mea- 
sured in units of Lj , such that there are d Lj lattice sites 
between the considered site and the step edge. The age 
r(d) of a position on the surface is proportional to dLj. 
The age will be measured in units of atoms added to 
the surface since the site has last been visited by a step. 
For N = 2 bunched steps moving together it amounts to 
r(d) = NLiLj d. 

The probability that there is no impurity at that site 
can be calculated as follows: The probability to turn a 
site into an impurity when a single atom is added to 
the surface is p = pjLiLj. Hence, the probability to be 
changed to an impurity after r(d) atoms have been added 
to the surface is 



Vi{d) = l-(l-p) 
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The latter approximation applies provided that 
LiLj/p ^S> 1. For the simulations shown in Fig. [2] 
one has LiLj/p — 25 000 such that this approximation is 
well justified. Moreover, in this situation the argument 
of the exponential function is small, such that the 
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FIG. 2: Impurity distributions (left) and color coded surface 
height profiles (right) at three different growth stages for a 
system with N = 2 steps which show step-pairing. The sys- 
tem parameters are p = 0.1; Li = 25; Lj = 100; J = 40; 
Ji m = 4; v — 2 x 10 6 . In the left panels red indicates impuri- 
ties and blue normal atoms in the topmost layer of the SOS 
representation of the lattice. The respective heights are given 
in the right panels with a color coding indicated by the bars 
on the far right. 

The upper panels show the initial configuration with two 
equidistant straight steps, and no impurities on the surface. 
The middle panels show the situation after the deposition 
of 0.2 ML. The impurities are roughly uniformly distributed, 
and the distances between the steps does not yet deviate much 
from the initial configuration. 

The lower panel shows the surface morphology after deposi- 
tion of 20 ML, where the steps sit right next to each other. 
Note that impurities at the terrace are not uniformly dis- 
tributed any more. There are more impurities ahead of the 
steps than behind the steps. 



distribution of impurities is approximately linear, 
Vi{d) = Npd, for Npd<^ 1. 



(8) 



Fig. [3] demonstrates that this prediction holds to a very 
good approximation for the time average over 8 ML. Plot- 
ting the ratio of the numerical values and the theoretical 
prediction (lower panel) shows that the agreement is bet- 
ter than 5% along the full width of the terrace. 

In contrast to time averages the distribution is very 
noisy, however, for snapshots of the distribution of impu- 
rities as shown by three different examples in the lower 
panel of Fig.[3J Indeed, for the considered deposition pro- 
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FIG. 3: Time average for the number of impurities as a func- 
tion of the distance d from the step for the simulation shown 
in Fig. [2] The upper panel shows a time average over a time 
span needed to grow 8 ML with numerical values indicated 
by red crosses and the theoretical prediction J7J by a dashed 
green line. The lower panel shows the ratio of the numerical 
values and the prediction (red crosses) together with those 
calculated from three different snapshots. The dashed green 
lines give an estimate of the width of the distribution for in- 
dividual snapshots (thick line) and a narrower band with a 
width that is smaller by a factor of five. 

cess the number of impurities in a row d is distributed ac- 
cording to a binomial distribution such that one expects 
to find Md — LiPiid) impurities in a row at a normal- 
ized distance d from the steps. The standard deviation 
should be LiVi(d) (1 — Vi(d)) such that the relative er- 
ror, which is indicated by a thick dashed green line in the 
lower panel of Fig. [3J takes the value 

SATd = ( 1-Vj(d) \ 1/2 
Nd V UViid) j ■ 

For small d hardly any atoms impinged on that part of 
the surface such that the probability V% is small. In that 
case the relative error exceeds 100%. From the perspec- 
tive of understanding the transport problem this is not 
so problematic, however, because the impurities hardly 
influence the diffusion on the surface in that case. It is 
much more important to observe that the relative error 
remains fairly large even far away from the steps, where 
the coverage of impurities is large. In the case of Fig. [3J 



where we consider expectation values for Li = 25 and the 
final coverage with impurities is close to 25% the relative 
error still amounts to (0.8/(25/5)) 1 / 2 = V0A6 = 0.4, 
and even for a coverage of 50% it only drops down to 0.2. 
Fluctuations of this magnitude are troublesome, when 
trying to describe the transport by a continuum model 
which does not take into account fluctuations of the dis- 
tribution of impurities. The problem can not be resolved 
by considering averages of larger Lj because the contin- 
uum description has to be based on local averages, and 
a one-dimensional description will only apply when the 
two-dimensional equations obtained in this manner are 
invariant under translation parallel to the steps. It is 
this latter property, however, which is lost when the fluc- 
tuations in the density of impurities are noticeable. 

The observation that the variance in the number of 
impurities in small neighborhoods of the lattice is large 
probably applies in general. This poses a major challenge 
to continuum models of step bunching where the presence 
of impurities is only taken into account as a modification 
of the diffusion coefficient, which itself depends on the 
expectation value for the number of impurities. In view 
of the large fluctuations of the distribution this might 
very well be a non-admissible oversimplification, which 
deserves a close inspection by comparison to numerical 
results in the following section. 

IV. IMPURITY-INDUCED STEP DYNAMICS 

A. Diffusion bias 

The key ingredient of the continuum theory developed 
in [9j is the dependence of the effective adatom diffusion 
coefficient D{&) on the local impurity coverage 9. As 
we are concerned here with impurity concentrations of 
6 ~ 0.1 or less, the leading term in an expansion in 9 is 
expected to suffice. We therefore write 

D{6) ~ D(0) (1 - a6). (9) 

For completely blocking barriers (Jj m — > oo) the coef- 
ficient is given bya = 7r — 1~2.14 [ijj]. In effective 
medium approximation [ill ] one obtains the simple ex- 
pression 

For J; m — > oo it yields a = 2, which is close to the exact 
result. 

Using the general formulae derived in Q, the adatom 
currents j± to the ascending and descending 
steps bordering a terrace of width w can be computed 
from Eq. ©. For perfectly absorbing steps one finds 
that j± — <&wp± , where the attachment probabilities p± 
are independent of the terrace width w, and given by 

^ ap ln(l — ap) ^ + ^ ^ 
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FIG. 4: Time evolution of the normalized step distance d 
for TV = 2 steps and different choices of Ji m : (a) attractive, 
Ji m = 4, (b) marginally stable, J lm — 1, and (c) repulsive, 
Jim = 0.2, interaction between the steps. The panels on the 
left show the evolution of an initial condition where the two 
steps are right next to each other, and the right panels that of 
an initial condition with equidistant steps. The other system 
parameters were fixed to the values L; = 25; Lj/N — 50; 
J = 40; v = 2 x 10 6 ; p = 0.1. 



p_ of the flux incident on the (trailing) terrace behind 
the step, of width x k — Xk-i- Thus we have 

^ = ~(1 - b)(x k+1 - x k ) + ~(1 + b)(x k - x k -x) (14) 

with the additional constraint that x k > x k ~i at all 
times. For convenience we use here dimensionless units 
where time is measured in units of the time scale $ _1 
needed to deposit a monolayer of new material, and 
length still in units of Lj. 



B. Step-step interactions 

Before turning to the discussion of our simulation re- 
sults, we need to address the role of repulsive step-step 
interactions that are usually added to the right hand 
side of (fT4"j) [l2l Hc| . As no direct step-step interactions 
are included in our KMC model, we only consider the 
well-known entropic interactions induced by collisions be- 
tween neighboring steps, which in turn are a consequence 
of thermal step meandering. Following [l7| . we estimate 
the distance between two such collisions along the trans- 
verse (i-) direction to be of the order of 
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where 8 is the step stiffness and w denotes the distance 
between the two steps. Clearly step collisions are irrele- 
vant as long as L c is larger than the lattice size Li parallel 
to the steps. We therefore conclude that the range of the 
repulsive step-step interactions in our simulations is lim- 
ited to step distances smaller than w c ~ {LiksT '/5) 1 / 2 . 
Using the expression 0, [l7| 



k R T 
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where the relation p + = 1— p_ ensures flux conservation. 
For small ap the attachment probabilities p± amount to 
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such that the effect of the impurities can be quantified 
by the diffusion bias parameter 



ap 

b = p--p+~ —. 
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The sign is chosen such that step bunching results for 

b>o mm. 

The following considerations will be based on a sim- 
ple one-dimensional, deterministic dynamical model for 
the step positions x k , k — 1,-,N measured along the 
j-direction. The speed of the A; th step is the sum of the 
fraction p + of the flux incident on the (leading) terrace 
in front of the step, of width Xk+\ — x k , and the fraction 



for the step stiffness in the SOS model at low tempera- 
tures, we find that w c ~ 2.8 for Li — 25 and J = 40. 
Thus the step-step interactions in our simulations are a 
purely local effect which merely ensures that steps cannot 
overtake each other. In this sense the situation is similar 
to that considered in (l9| within a one-dimensional model 
with hard core step-step interactions. 



C. Stability of step pairs 

Consider first a system of two steps, with normalized 
step distances d(t) and 1 — d(t) (Fig. 2]). It follows from 
(fl"4|) that d(t) evolves according to 



d = b(2d- 1). 



(17) 



The equidistant fixed point d = 1/2 is unstable (stable) 
for b > (b < 0). For b > the steps collide (d — > 
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FIG. 6: The probability densities to find two adjacent steps 
at a distance w for systems with TV = 3 (top) and TV = 8 
(bottom) steps, respectively. For TV = 3 the lines (a)-(c) 
refer to the data shown in the respective panels of Fig. [5] and 
the simulation with TV = 8 was run with Li = 20; Lj/N = 50; 
J = 40; v = 8 x 10 6 ; p = 0.1; and (a) J im = 4, (b) J im = 1, 
(c) Ji m = 0.2. The slightly smaller width Li = 20 and larger 
v was chosen to minimize the impact of island formation. As 
also suggested by the histograms this change has no significant 
impact on the step interaction. 



FIG. 5: Time evolution of normalized step distances cfe and 
dz for TV = 3 steps and different choices of Ji m : (a) attrac- 
tive, Jim = 4, (b) marginally stable, Ji m = 1, and (c) re- 
pulsive, Jim = 0.1 interaction between the steps. The other 
parameters are Li = 25; Lj/N = 50; J = 40; v = 3 x 10 6 ; 
p — 0.1. When impurities induce repulsion (c), the steps 
remain well-separated, just as in the simulations with only 
two steps. However, for three steps a smaller value of Ji m 
(Jim = 0.1 rather than J lm = 0.2 displayed in Fig. [4] for 
two steps) was needed to clearly show this effect. In the 
marginally stable case (b) the distances fluctuate showing reg- 
ularly looking oscillations. They arise due to coupling of step 
velocities because of shared neighboring terraces. In this case 
the steps approach each other much closer than in the repul- 
sive case. However, still they always remain well-separated 
due to entropic repulsion, which is always present. Finally, in 
the case of impurity-induced attraction between the steps (a), 
one clearly sees bunches of two steps, which separate however 
when the third step approaches. 



0) in finite time, while for b < a pair of close steps 
separates and approaches the fixed point exponentially 
at rate 2b. Using the expression (|13[) with a given by 
Eq. (JTUJ) the half-time of the decay into the stable states, 



tj/2 = In 2/26 is about 16 ML. This value is consistent 
with but somewhat larger than the time scale observed 
in the simulations Fig. 0|a) and (c) . A possible source of 
this deviation is the fact that the impurity concentration 
profile may not have reached stationarity on the time 
scale of step motion. 



D. Stability of small bunches 

Consider next a system of three steps, where d,2 and 
^3 denote the normalized distance from the first to the 
second and the first to the third step, respectively. Fig. [5] 
shows the evolution of c?2 and dj, as a function of time t 
measured in units of deposited ML. 
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In the absence of impurity-induced step interactions 
(Jim = 1; Fig. [5(b)) and for repulsive impurity-induced 
step interactions (</;,„ = 0.2, Fig. 0c)) the system be- 
haves very similar to the one with only two steps. In- 
deed, histograms for the probability to find a certain dis- 
tance w between adjacent steps in Fig. [6] consistently 
show sharply peaked distributions around the the aver- 
age terrace size w = 50, while the distribution is broad 
with a maximum at this value in the neutral case. As ex- 
pected for a system with attractive interactions between 
the steps, the distribution P(w) has considerably more 
weight for small w. However, surprisingly, it does not 
decay but saturates at fairly large constant background 
extending till w = 100. 

The origin of this background becomes clear from 
inspection of the time traces of c?2 and ds shown in 
Fig. 02a) . The simulation shows the transient formation 
of step pairs which exchange partners at regular intervals. 
However, no stable step triplets are formed. To see how 
this follows from the dynamical equations (fT4|) , suppose 
a triplet of three nearby steps has formed, such that the 
step positions satisfy xi — x%, X3 — xi -C 1, i.e., they are 
both much smaller than the system size Lj. Then step 3 
has a large terrace of size ~ 1 in front of it, and it moves 
at speed (1 — b)/2. Step 2 is surrounded by small terraces 
and moves very slowly, and step 1 is constrained by the 
no-passing condition to move at the same speed. Thus 
for every b < 1 step 3 will detach from the triplet. 

We conclude that, for any b < 1, step triplets and 
larger bunches are unstable against the emission of single 
steps. Bound states of steps moving together at constant 
speed [20] can form when b > 1 [T4], but this condition 
obviously cannot be reached in the growth model con- 
sidered here. A detailed analysis of the equations (fl"4")) 
shows that step emission will continue until the number 
of free steps between two bunches (or, equivalently, be- 
tween one bunch and its periodic image) has reached the 
steady state value [l2| 

N f ~±lnN, (18) 

where N denotes the total number of steps in the bunch 
and on the terrace. As Nf cannot exceed N, we conclude 
that a stable bunch can form only if the number of steps 
satisfies the condition 

3bN/lnN > 1. (19) 

With the value of b ~ 0.02 obtained for J lm = 4 and 
p = 0.1 this implies N > 65, which (given the constraints 
on the systems parameters described above) exceeds our 
computational capacities. Indeed, simulations conducted 
with systems containing up to 20 steps confirm that step 
bunching remains a transient phenomenon, even when a 
bunched initial configuration is chosen (Fig. [7]) . 
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FIG. 7: Evolution of a bunch of eight (top) and twenty (bot- 
tom) steps in space-time representation. The simulation with 
N — 8 was run with the same parameters as in Fig. [6] (bot- 
tom^) except for a larger value v = 2.4 x 10 7 , and the one 
with N = 20 with Li = 20; Lj /N = 50; J = 40; J im = 4; 
v = 3 x 10 7 and p = 0.1. 



E. Stability of large bunches 

The agreement between the predictions and the results 
of the kinetic Monte-Carlo simulations indicates that the 
stability of bunches can be captured based on the one- 
dimensional growth model. To illustrate the formation 
of stable bunches with increasing N we therefore relied 
on simulations of a one-dimensional stochastic growth 
model described in detail in [l9j]. In this model parti- 
cles are deposited onto a one-dimensional vicinal surface 
and transferred instantaneously (without explicit diffu- 
sion) to the ascending or descending step with probabil- 
ities p + and p_ , respectively. Since steps are allowed to 
coalesce (but not to pass each other), the formation of 
bunches can easily be followed by monitoring the maxi- 
mal step height /i max (the largest nearest neighbor height 
difference) in the system. In Fig. [5] we show results ob- 
tained for b = 0.02 on a lattice of Lj = 1000 sites. For 
iV = 20, the situation corresponding to Fig. [7(b), the 
maximal step height remains below 2, showing that only 
step pairs exist, whereas for N — 100 a stable bunch 
forms that contains almost half of the steps in the sys- 
tem. 
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FIG. 8: Simulations of step bunching in a one-dimensional 
model. The top and middle panels show the time evolution 
of the maximal step height for systems containing N = 20 
and N— 100 steps, respectively. The respective data points 
in these plots correspond to initial conditions with a single 
large step (red) and an equidistant step train (green), which 
were averaged over 100 independent runs. Beware of the dif- 
ferent time scales of these two plots. The inset in the upper 
graph shows the ratio of the data of the ID model (red curve) 
to the maximal height of a bunch in the corresponding KMC 
simulation Fig. [7] The bunch size was calculated by deter- 
mining the maximal number of steps in a window of fixed 
width d = 0.03, d = 0.06, and d = 0.09 (from bottom to top) 
roughly corresponding to the range of the effective hard-core 
repulsion Nw c /Lj. The lowermost panel shows the deviation 
Ahj = hj — 10 5 of the height hj from the average height after 
deposition of 10 5 ML. It is the final configuration of one of the 
simulations run to generate the graph in the middle panel. 



V. CONCLUSIONS 

We summarize the main achievements of this work: 

(i) We have verified by KMC simulations that an 
impurity-induced increase of the transition energy be- 
tween neighboring sites, a purely kinetic effect, is a possi- 
ble source of step bunching. By focusing on the behavior 
of pairs of steps we have explicitly determined the sign 
and strength of the impurity-induced step interactions, as 
quantified by the asymmetry parameter b. The adatoms 
can make use of massive fluctuations of the spatial and 
temporal distribution of the impurities (Fig. [3|) in order 
to find optimal paths to the binding sites at the steps. As 
a consequence an estimate of the order of magnitude of b 
properly has to account for the diffusion of the adatoms 
in a 2d disordered arrangement of impurities, which is 
very different from the Id setting used in [9(. 

(ii) The step bunching observed in our work is substan- 
tially weaker than that found in simulations of the SiC 
system Q, where the impurities were assumed to affect 
the adatom binding energies. This probably indicates 
that the effective asymmetry b is larger for energetic im- 
purities. Unfortunately, although the theory of [§| cor- 
rectly predicts step bunching for impurities that lower 
the binding energy, the magnitude of the effect depends 
on the precise boundary conditions at the steps, which 
do not easily translate into the two-dimensional KMC 
setting. 

(iii) Despite the presence of impurity-induced attrac- 
tive step interactions, triplets and larger assemblies of 
steps are not necessarily stable when b is small: Instead of 
agglomerating into macroscopic step bunches, the steps 
display a peculiar dynamical pattern characterized by 
transient step pairs that exchange partners much like in 
a folk dance. Up to now this effect has gone unnoticed, 
because previous KMC simulations of step bunching dur- 
ing growth have generally considered situations where 
the effective attachment asymmetry b is of order unity 
[2IL I22L I23I ] . The lack of stability of the step assemblies 
was explained based on a recently developed determin- 
stic theory [12j. It allows us to predict that bunches can 
form only when the number of steps exceeds the bound 

(iv) For real surfaces there is no restriction on the total 
number of steps. Nevertheless, it is highly improbable to 
observe bunching in systems with small b. Our simula- 
tions for the one-dimensional model (Fig. [5]) show that for 
small b bunches evolve only after exceedingly long times 
even when the bound (fl9|) is satisfied: In order to see step 
bunching one has to wait for a fluctuation nucleating a 
bunch with a minimal size given by (|19p. 

For applications the most noticeable consequence of 
our study is that the rapid formation of large step 
bunches seen experimentally in the growth of SiC [8| 
cannot be explained only in terms of kinetic impurities. 
Some coupling to the adatom binding energy must also 
be involved. 
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